Análisis Multivariante
Profesor: Laura Durán
Alumnos: Steven Allus, Antonio López
En esta actividad, exploraremos el Análisis de Componentes Principales (PCA) aplicando técnicas estadísticas en R.
La tarea consiste en llevar a cabo un Análisis de Componentes Principales a partir de un estudio en el que se ha medido la expresión de los genes sobre distintos tumores. Los datos están almacenados en el archivo genes_tumores.csv que se puede descargar desde el aula virtual. Esta tabla contiene los valores de expresión de 12 genes medidos en 10 tumores distintos.
El objetivo del análisis será averiguar si de los 10 tumores analizados existen algunos con características comunes en cuanto a la expresión de sus genes. Los resultados del análisis deberán representarse en un biplot en el que las variables (genes) se representen como vectores y los individuos (tumores) como puntos.
Pueden seguirse los pasos explicados en clase u otros alternativos que conduzcan a las mismas conclusiones.
El trabajo es en parejas
Aspectos formales de las tareas
Hay que entregar un informe con el código en R utilizado, las explicaciones pertinentes para cada paso del análisis y las conclusiones derivadas del análisis. Se recomienda usar RMarkdown.
Los Componentes Principales son combinaciones lineales de las variables originales diseñadas para maximizar la varianza explicada.
Datos con los que trabajar: Diez tumores con datos de doce genes.
Tenemos que reducir las variables de los tumores para representarlas visualmente siempre perdiendo información.
A partir de la combinación lineal de las doce variables obtendremos una reducción de estas que nos permitirá representarlas gráficamente.
Estas variables representadas gráficamente se llaman Componentes Principales.
Vamos a representar un eje de coordenadas utilizando subconjuntos de datos de los genes. Cada conjunto incluye tres genes, lo que permite una visualización tridimensional intuitiva. Los datos están representados con un gráfico interactivo que muestra cómo se distribuyen las muestras de tumores en función de las combinaciones lineales de los genes. Esta representación facilita observar similitudes o diferencias entre tumores en relación con los genes seleccionados.
# Establecer el repositorio de CRAN para evitar errores
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Seleccionar directorio de trabajo
setwd("/Users/stevenallus/Downloads")
# Instalar y cargar paquetes necesarios
if (!require(boot)) install.packages("boot")
library(boot)
paquetes <- c("logisticPCA", "FactoMineR", "factoextra", "plotly", "readxl")
paquetes_faltantes <- paquetes[!paquetes %in% installed.packages()[,"Package"]]
if (length(paquetes_faltantes) > 0) install.packages(paquetes_faltantes)
# Cargar las librerías
lapply(paquetes, library, character.only = TRUE)
# Importar datos desde un archivo CSV
genes_tumores <- read.csv("genes_tumores.csv")
genes_tumores <- genes_tumores[, -1]
head(genes_tumores)
## gene1 gene2 gene3 gene4 gene5 gene6
## 1 0.9323190 1.2113278 0.5383686 -0.3804224 -0.5994117 -0.46955658
## 2 1.0613959 0.7762260 0.5098463 0.1574786 -0.4035551 -0.50984955
## 3 0.4637662 1.1793956 1.1355446 0.1548933 -0.7319309 -0.42653172
## 4 -1.6228992 -0.7777483 -0.7868679 -1.1771294 -0.3007739 0.11591201
## 5 -1.1035586 -0.9511971 -0.9131699 -1.0484974 -0.3974149 -0.35492610
## 6 -1.0950318 -0.9722049 -0.9504966 -0.8261171 -0.5278552 -0.06271155
## gene7 gene8 gene9 gene10 gene11 gene12
## 1 -1.06758821 -1.23610151 0.1145246 0.8975893 1.4624217 1.1259876
## 2 -1.31883371 -0.87248832 -0.2017863 1.5210825 1.1883528 0.9711843
## 3 -1.03071756 -1.05634140 -0.2747214 1.0444794 0.7132720 0.7550505
## 4 0.19523024 0.15654077 0.7884622 0.9194994 1.0616514 0.9178455
## 5 0.23230945 -0.22601482 1.0759271 1.1158295 0.8121759 1.1303828
## 6 -0.04847066 0.05124849 1.3769491 0.8643747 0.9175155 0.6019269
Este código escala todas las variables para conseguir las componentes principales y nos muestra un gráfico.
El escalado (normalización) asegura que todas las variables tengan la misma importancia en el análisis, independientemente de sus unidades o magnitudes. Esto se hace transformando cada variable para que tenga una media de 0 y una desviación estándar de 1.
Sin el escalado, las variables con valores más grandes (como pesos en kilogramos) dominarían a las variables con valores más pequeños (como proporciones o porcentajes).
Las componentes principales son combinaciones lineales de las variables originales, y el escalado garantiza que estas combinaciones no estén sesgadas por las escalas originales de las variables. Esto permite que el análisis se enfoque únicamente en las relaciones entre las variables y no en sus magnitudes.
# PCA -> datos escalados, con plot:
pca_genes <- PCA(genes_tumores, scale.unit = T, graph = T)
Gráfico de PCA mostrando la reducción de dimensiones.
Conclusiones:
Esto significa que los dos primeros componentes principales explican el 96.45% de la información contenida en los datos originales.
Ordena las dimensiones de mayor a menor valor de eigenvalue.
Eigenvalues o valores propios: mide la cantidad de varianza explicada por esa componente.
# str(pca_coronaria) para ver la estructura del objeto
summary(pca_genes)
##
## Call:
## PCA(X = genes_tumores, scale.unit = T, graph = T)
## Eigenvalues
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 Dim.7
## Variance 8.183 3.391 0.132 0.095 0.070 0.053 0.050
## % of var. 68.191 28.260 1.098 0.789 0.579 0.443 0.417
## Cumulative % of var. 68.191 96.451 97.549 98.338 98.918 99.361 99.778
## Dim.8 Dim.9
## Variance 0.017 0.009
## % of var. 0.144 0.079
## Cumulative % of var. 99.921 100.000
##
## Individuals
## Dist Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr cos2
## 1 | 2.713 | -1.819 4.043 0.450 | 1.903 10.682 0.492 | -0.141 1.513 0.003 |
## 2 | 2.594 | -1.543 2.910 0.354 | 1.928 10.962 0.552 | 0.485 17.858 0.035 |
## 3 | 2.494 | -1.071 1.401 0.184 | 2.085 12.822 0.699 | -0.654 32.420 0.069 |
## 4 | 3.363 | -2.178 5.798 0.419 | -2.485 18.207 0.546 | -0.220 3.678 0.004 |
## 5 | 3.464 | -2.656 8.622 0.588 | -2.160 13.763 0.389 | 0.186 2.614 0.003 |
## 6 | 3.391 | -2.363 6.824 0.486 | -2.352 16.316 0.481 | 0.096 0.694 0.001 |
## 7 | 4.218 | 4.100 20.541 0.945 | -0.717 1.518 0.029 | -0.470 16.784 0.012 |
## 8 | 4.199 | 4.168 21.230 0.985 | -0.271 0.216 0.004 | -0.032 0.079 0.000 |
## 9 | 4.711 | 4.662 26.566 0.980 | -0.215 0.136 0.002 | 0.515 20.123 0.012 |
## 10 | 2.697 | -1.300 2.065 0.232 | 2.284 15.377 0.717 | 0.236 4.238 0.008 |
## Variables (the 10 first)
## Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr cos2
## gene1 | 0.530 3.434 0.281 | 0.808 19.260 0.653 | 0.198 29.684 0.039 |
## gene2 | 0.577 4.066 0.333 | 0.800 18.876 0.640 | -0.093 6.571 0.009 |
## gene3 | 0.595 4.321 0.354 | 0.787 18.268 0.620 | -0.100 7.620 0.010 |
## gene4 | 0.910 10.117 0.828 | 0.372 4.079 0.138 | 0.019 0.279 0.000 |
## gene5 | 0.949 11.014 0.901 | -0.220 1.433 0.049 | 0.191 27.663 0.036 |
## gene6 | 0.917 10.278 0.841 | -0.378 4.217 0.143 | -0.040 1.211 0.002 |
## gene7 | 0.707 6.114 0.500 | -0.674 13.390 0.454 | -0.102 7.833 0.010 |
## gene8 | 0.781 7.449 0.610 | -0.611 11.002 0.373 | 0.092 6.401 0.008 |
## gene9 | -0.867 9.185 0.752 | -0.449 5.952 0.202 | 0.051 1.980 0.003 |
## gene10 | -0.963 11.328 0.927 | 0.198 1.154 0.039 | 0.064 3.090 0.004 |
Los eigenvalues miden la varianza explicada por cada componente.
plot(
pca_genes$ind$coord[, 1], # Coordenadas en Dim1
pca_genes$ind$coord[, 2], # Coordenadas en Dim2
xlab = paste0("Dim 1 (", round(pca_genes$eig[1, 2], 2), "% varianza)"),
ylab = paste0("Dim 2 (", round(pca_genes$eig[2, 2], 2), "% varianza)"),
main = "Proyección en Dim1 y Dim2", # Título del gráfico
pch = 19, # Tipo de punto
col = "blue", # Color de los puntos
xlim = c(min(pca_genes$ind$coord[, 1]) - 1, max(pca_genes$ind$coord[, 1]) + 1),
ylim = c(min(pca_genes$ind$coord[, 2]) - 1, max(pca_genes$ind$coord[, 2]) + 1) )
# Eje coordenadas
abline(h = 0, v = 0, col = "black", lty = 2)
# Etiquetas
text(
pca_genes$ind$coord[, 1],
pca_genes$ind$coord[, 2],
labels = rownames(pca_genes$ind$coord),
pos = 4, # Posición del texto
cex = 0.7 # Tamaño del texto
)
Resultados —
Proyección de los datos en Dim1 y Dim2.
# Varianza explicada por cada dimensión:
fviz_eig(pca_genes, addlabels = TRUE, ylim = c(0, 100))
library(ggplot2)
library(plotly)
# Crear el data frame
varianza <- data.frame(
Dimensión = 1:length(pca_genes$eig[, 1]),
Varianza = pca_genes$eig[, 2],
Acumulada = cumsum(pca_genes$eig[, 2])
)
# Gráfico con ggplot2
graf_varianza <- ggplot(varianza, aes(x = Dimensión)) +
geom_bar(aes(y = Varianza), stat = "identity", fill = "steelblue") +
geom_line(aes(y = Acumulada), color = "red", linewidth = 1) +
geom_point(aes(y = Acumulada), color = "red", size = 2) +
labs(
title = "Varianza Explicada y Acumulada",
x = "Dimensión",
y = "Porcentaje de Varianza"
) +
theme_minimal()
# Exportar el gráfico a plotly
ggplotly(graf_varianza)
# Librerías necesarias
library(FactoMineR)
library(factoextra)
library(plotly)
# Crear el biplot con FactoMineR
biplot_pca <- fviz_pca_biplot(
pca_genes,
repel = TRUE, # Evitar superposición de etiquetas
col.var = "blue", # Color de las flechas de variables
col.ind = "black", # Color de los puntos de individuos
arrowsize = 1, # Tamaño de las flechas
pointsize = 3, # Tamaño de los puntos
geom = c("point", "text"),
title = "Biplot: Genes (Flechas) y Tumores (Puntos)"
)
# Exportar el biplot a plotly para interactividad
ggplotly(biplot_pca)
Vamos a validar la estabilidad y consistencia del análisis realizado. Para ello haremos empleo la técnica de remuestreo (bootstrap) para evaluar la estabilidad de los componentes principales. Aquí está la propuesta:
# Función para remuestrear los datos y recalcular el PCA
pca_bootstrap <- function(data, indices) {
# Subconjunto de datos remuestreados
data_boot <- data[indices, ]
# Realizar PCA con los datos remuestreados
pca_result <- PCA(data_boot, scale.unit = TRUE, graph = FALSE)
# Concatenar las coordenadas de PC1 y PC2 en un vector
c(pca_result$ind$coord[, 1], pca_result$ind$coord[, 2])
}
# Configurar remuestreo bootstrap
set.seed(123) # Para reproducibilidad
bootstrap_results <- boot(
data = genes_tumores, # Datos originales
statistic = pca_bootstrap, # Función definida arriba
R = 1000 # Número de remuestreos
)
# Resumen de resultados
print(bootstrap_results)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = genes_tumores, statistic = pca_bootstrap, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -1.8188739 2.0286898 2.995240
## t2* -1.5431329 1.6136273 2.914554
## t3* -1.0708431 0.8549935 2.761629
## t4* -2.1782566 2.0206275 2.808105
## t5* -2.6561326 2.7046518 2.950829
## t6* -2.3630028 2.3056353 2.877474
## t7* 4.0998488 -4.1402529 2.853162
## t8* 4.1679928 -4.0180110 2.850127
## t9* 4.6624718 -4.7633048 2.823456
## t10* -1.3000715 1.3933435 2.936113
## t11* 1.9032783 -1.8512380 1.749042
## t12* 1.9280936 -1.9604068 1.804571
## t13* 2.0852384 -2.1047527 1.862594
## t14* -2.4848021 2.5909666 1.831133
## t15* -2.1603981 2.1103791 1.762714
## t16* -2.3522557 2.3860561 1.771807
## t17* -0.7174422 0.7045590 1.848779
## t18* -0.2705048 0.2949778 1.787310
## t19* -0.2147895 0.1865351 1.854436
## t20* 2.2835821 -2.3570762 1.856363
# Visualización de la variación en las componentes principales
plot(bootstrap_results, index = 1, main = "Variación en PC1 (Bootstrap)")
plot(bootstrap_results, index = 2, main = "Variación en PC2 (Bootstrap)")
De los resultados del análisis de bootstrap para los tumores (t1* a t19*), se pueden extraer las siguientes conclusiones basadas en los valores presentados y los gráficos (histograma y QQ-plot):
El análisis bootstrap confirma que: 1. Los patrones de agrupación observados en el biplot son consistentes, especialmente para tumores con menor bias y std. error. 2. Algunos tumores atípicos (como t7, t8, t9* y t14*), aunque claramente identificados en el biplot, son más sensibles al remuestreo, lo que sugiere que sus posiciones deben interpretarse con precaución. 3. Tumores centrales en el biplot tienen posiciones robustas y estables, reforzando su relevancia en el análisis.